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For systems of interacting, ultracold spin-zero neutral bosonic atoms, harmonically trapped and 
subject to an optical lattice potential, we derive an Extended Bose Hubbard (EBH) model by devel- 
oping a systematic expansion for the Hamiltonian of the system in powers of the lattice parameters 
and of a scale parameter, the lattice attenuation factor. We identify the dominant terms that need 
to be retained in realistic experimental conditions, up to nearest-neighbor interactions and nearest- 
neighbor hoppings conditioned by the on site occupation numbers. In mean field approximation, 
we determine the free energy of the system and study the phase diagram both at zero and at finite 
temperature. At variance with the standard on site Bose Hubbard model, the zero temperature 
phase diagram of the EBH model possesses a dual structure in the Mott insulating regime. Namely, 
for specific ranges of the lattice parameters, a density wave phase characterizes the system at integer 
fillings, with domains of alternating mean occupation numbers that are the atomic counterparts of 
the domains of staggered magnetizations in an antiferromagnetic phase. We show as well that in 
the EBH model, a zero-temperature quantum phase transition to pair superfluidity is in princi- 
ple possible, but completely suppressed at lowest order in the lattice attenuation factor. Finally, 
we determine the possible occurrence of the different phases as a function of the experimentally 
controllable lattice parameters. 



I. INTRODUCTION 

In the last years we have witnessed a spectacular ac- 
celeration in the experimental manipulation of neutral 
atoms in optical lattices 0, 0, IS that has opened the 
way to the simulation of quantum complex systems of 
condensed matter physics, such as high-X^ superconduc- 
tors, Hall systems, and superfluid 4 _ffe, thanks to the ex- 
treme flexibility and controllability of such atom-optical 
systems [Bj. Optical lattices are stable periodic arrays 
of microscopic potentials created by the interference pat- 
terns of intersecting laser beams 6]. Atoms can then 
be confined in different lattice sites, and by varying the 
strength of the periodic potential it is possible to tune 
the interatomic interactions with great precision and to 
enhance them well into the regime of strong correlations, 
even in the dilute limit. The transition to a strong cou- 
pling regime can be realized by increasing the depth of 
the lattice potential wells, a quantity that is directly pro- 
portional to the intensity of the laser light which, in turn, 
is an experimental parameter that can be controlled with 
great accuracy. For this reason, besides the fundamental 
interest for the investigation of quantum phase transi- 
tions I d |c |l9 |.|l0| and other collective quantum phenom- 
ena Il2l Il3l If"H [TEL [Lsl ] , optical lattices have become 
an important tool in applications ranging from laser cool- 
ing hjj to quantum control and information processing 
[18l Il9| , and quantum computation |20l l2lL |22j, I23L 12 1| . 

The theory of neutral bosonic atoms in optical lattices 
has been originally developed, in its simplest framework, 
by assuming that the atoms are confined to the lowest 
Bloch band of the periodic optical potential Hj||26]. It is 
then straightforward to show that in this approximation 



the system is very well described by the standard, on site 
Bose Hubbard model. In such a model a superfluid-Mott 
insulator (SF-MI) transition is predicted to occur when 
the energy gap between the local ground state and the 
first excited levels becomes comparable to the hopping 
energy between adjacent lattice sites 0, IE IH E3 • More- 
over, no major qualitative changes appear when higher 
excited energy levels of the external trapping potentials 
are considered |27| . This prototypical quantum phase 
transition can be realized experimentally, for instance by 
manipulating the strength of the lattice potential, which 
results in a controlled change of the kinetic (hopping) 
energy term Following this kind of approach, the 
superfluid-Mott insulator quantum phase transition has 
been realized in a series of beautiful experiments, by load- 
ing an ultracold atomic Bose-Einstein condensate in a 
three-dimensional optical lattice |4j. At finite temper- 
ature and in the strong coupling regime, the quantum 
phase transition is smeared in a classical transition from 
a completely disordered regime to an ordered, superfluid 
phase (SF) [lot l27| . Finally, if more complicated geomet- 
rical settings such as optical superlattices, networks, and 
graphs, are considered, interesting unconventional inter- 
plays between the different quantum phases may appear 

Mm in mm. 

In the on site Bose Hubbard (BH) model all the terms 
deriving from two-body interactions between the bosonic 
atoms are neglected except the local ones taking place 
on a same site of the optical lattice. This is often an 
excellent approximation. However, microscopic interac- 
tions are in general of finite range, and it may be in- 
teresting to address the problem of the physical picture 
that emerges if one takes into account two-body inter- 
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actions between sites at non- vanishing distance and/or 
hopping amplitudes beyond nearest-neighbor sites. This 
problem has been extensively studied in the case of 
charged fermions, which are endowed with very long 
range Coulomb interactions, usually by introducing ex- 
tended Fermi Hubbard models with nearest-neighbor in- 
teraction terms [HIM SIM El- Extended Fermi Hub- 
bard models are prototypical in theoretical condensed 
matter physics, because they are believed to be more re- 
alistic approximations to the true inter-electronic inter- 
actions and moreover exhibit richer phase diagrams when 
compared to the standard Hubbard model, especially in 
low spatial dimensions j^l, |M 0] . 

In the present paper, we extend this kind of analysis to 
the bosonic case, and we introduce and study a type of 
extended Bose Hubbard model for the description of in- 
teracting bosons in regular lattices. Our study, although 
general, will be mainly concerned with the analysis of 
realistic physical systems such as dilute ensembles of in- 
teracting spin-zero bosonic neutral atoms subject to an 
optical lattice potential and spatially confined by a slowly 
varying harmonic trapping. Confining our attention to 
the lowest band of the optical lattice potential, we con- 
struct an effective extended Bosc Hubbard Hamiltonian, 
that emerges when we take into account the terms deriv- 
ing from the boson-boson interactions both on the same 
lattice site and on pairs of nearest neighbors. Our aim is 
then to analyze the phase transitions that may occur in 
such systems both at zero and at finite temperature, and 
determine the general features of the associated phase 
diagrams. 

The paper is organized as follows: In Sec. II we set the 
general notations, the needed tools, and derive the model 
Hamiltonian from the underlying microscopic dynamics 
in second quantization. We identify an exponential lat- 
tice scale factor, the lattice attenuation factor that gives 
rise to a natural expansion parameter for the Hamilto- 
nian of the system. In this way we are able to determine 
all the dominant terms that need to be retained in a 
power series expansion of the different energy contribu- 
tions, that add to the terms of the standard Bose Hub- 
bard (BH) model Hamiltonian. These additional terms 
can be easily identified and ordered in increasing pow- 
ers of the lattice attenuation factor and include nearest- 
neighbor interactions, nearest-neighbor hoppings of sin- 
gle atoms conditioned by the on site occupation number, 
and nearest-neighbor hoppings of atomic pairs (bosonic 
pairs). The ensuing Hamiltonian defines a type of ex- 
tended Bose Hubbard (EBH) model of interacting bosons 
in the lowest Bloch band. Other types of EBH mod- 
els can be obtained by considering interactions with far 
neighbors, and can be related to the fractional quantum 
Hall effect |4l|. Other different types of Hubbard-like 
models can also be obtained, at least for one-dimensional 
fermionic systems, by implementing appropriate mode 
expansions of the second-quantized fermionic field oper- 
ators, as recently shown by Massel and Penna ^1 ■ 

The quantum phase diagrams of different types of EBH 



models have been extensively investigated by means of 
quantum Monte Carlo techniques. In these studies the 
possible existence of density wave and supersolid phases 
in 1 and in 2 spatial dimensions has been carefully dis- 
cussed 0, Bl- hi t ne present work we set up a for- 
malism that allows to study in detail the explicit depen- 
dence of EBH models on the physical parameters of re- 
alistic systems of ultracold neutral atoms in optical lat- 
tice potentials and external confining harmonic poten- 
tials. Equipped with these tools, we can then establish 
the range of values of the Hamiltonian parameters that 
should be reached experimentally for the possible obser- 
vation of the new quantum phases predicted by EBH 
models, such as the density wave Mott insulating regime. 
Our treatment does not involve issues of metastability 
and lifetime for states defined on excited energy bands; it 
allows the study of the phase diagram both at zero and at 
finite temperature; and, finally, is based on a systematic 
expansion that connects any EBH model to the underly- 
ing microscopic many-body dynamics in second quantiza- 
tion. Moreover, it can be in principle extended to more 
general situations, including multi-species bosonic sys- 
tems and Fermi-Bose mixtures of bosonic and fermionic 
atoms. Other generalizations of the on site Hubbard 
model are of course possible by including Bloch bands 
of higher order in the description. Two such possible 
EBH models with inter-band hoppings and interactions 
have been recently investigated in connection with pro- 
posed schemes to generate metastable excited states of 
cold atoms in optical lattices [45|, |4(j . 

In Sec. Ill, we analyze the phase diagram of the ex- 
tended Bose Hubbard model (EBH) in the standard 
grand canonical formalism. We first consider the struc- 
ture of the quantum phases at zero temperature, when 
all the kinetic terms vanish (strong coupling regime) . Re- 
markably, in this case the model is mapped into a quan- 
tum antiferromagnetic Ising model in the presence of an 
external magnetic field. As it is well known, this model 
undergoes a quantum transition between a ferromagnetic 
and an antiferromagnetic phase |47^ . As a consequence, 
in a small range of parameters (very strong coupling 
regime) for a system of bosons in an optical lattice, there 
exists a new insulator phase in which the atomic den- 
sity is not constant on each lattice site. This new quan- 
tum phase can be seen as a Density Wave Mott Insulator 
(DWMI) , and is of course absent in the on site BH model 
that can sustain only a Pure Mott Insulator (PMI) phase 
at constant density. In the DWMI phase, there appear 
two alternating domains (sublattices) of different mean 
on site occupation numbers, that are the atomic coun- 
terparts of the domains of different staggered magnetiza- 
tion in the antiferromagnetic phase of the quantum Ising 
model. Next, we reintroduce the kinetic terms in the 
EBH Hamiltonian by resorting to Bogoliubov-like and 
mean field approximations. In this framework three real- 
valued order parameters emerge: the conventional single 
boson (single atom) superfluid order parameter, a new 
bosonic pair superfluid order parameter, and finally the 
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mean number of bosons (atoms) per site. This parame- 
ters are not all independent from each other, due to the 
existence of physical constraints on thermodynamic sta- 
bility and on the PMI-DWMI phase separation. 

In Sec. IV, we present qualitative and analytical stud- 
ies of the different possible phases in the presence of non 
vanishing kinetic energy terms. In this case, minimiza- 
tion of the free energy with respect to the different or- 
der parameters yields that the only superfluid transition 
which can occur is the one ruled by the single atom SF 
order parameter. We analyze the behavior or the SF- 
PMI and the SF-DWMI quantum phase transitions at 
zero temperature, and the behavior of the SF single bo- 
son order parameter at finite temperature. We study the 
critical temperature of the disordered-SF phase transi- 
tion as a function of the filling factor and of the lattice 
depth. Then, starting from finite temperatures, we an- 
alyze the possibility to recover the SF-PMI and the SF- 
DWMI quantum phase transitions by determining the 
behavior of the critical energy gap in the limit of vanish- 
ing temperature, as a function of the lattice depth and 
for different values of the trapping frequencies. By com- 
paring the phase diagrams obtained for different values 
of the Hamiltonian parameters, we discuss the possibility 
of observing experimentally the zero-temperature transi- 
tion to the SF-DWMI phase in systems of neutral atoms 
loaded in optical lattice potentials. Finally, in Sec.V we 
give some concluding comments on the obtained results 
and discuss some possible directions for future research. 



II. GENERAL SETTING 

The microscopic Hamiltonian for an ensemble of 
bosonic atoms that are confined by a slowly varying ex- 
ternal harmonic trapping potential and subject to an ad- 
ditional optical lattice can be written as 

H = f + V + W , (1) 

where T is the kinetic energy term that reads 



f = / dr¥(r)V 2 i>{r) 

2m J 



(2) 



with ^(r) being the bosonic annihilation field operator at 
point r. On the other hand V represents the contribution 
of the external potential to the energy: 

V = J dr¥(f){V H {r f ) + V opt (f))i>(f) . (3) 

In concrete situations, spatial confinement of the atoms 
is provided by the quadrupolar anisotropic trapping mag- 
netic field that leads to an harmonic trapping potential 
of the form 



where lo is the frequency associated to the harmonic trap 
in the x direction and A is the anisotropic coefficient, 
i.e. the ratio between the frequency in the yz plane 
and the frequency in the x direction. The most com- 
mon experimental settings are realized in the so-called 
"cigar-shaped" configuration (A ^> 1). Correspondingly, 
the second contribution to the potential energy is a 1-D 
periodic potential needed to realize the optical lattice 
along the axis of the "cigar" : 



V opt (x) = Vosin 2 (^) 



(5) 



where Vb is the maximum amplitude of the light shift 
associated to the intensity of the laser beam and a is 
the lattice spacing related to the wave vector k of the 
standing laser light by k — ir/a. 

Finally, the third contribution is the local, contact two- 
body interaction 



W 



9bb 



(6) 



in which the interaction coupling gsB = 4-7r?i obb/jti 
where m is the atomic mass and ass is the atom-atom 
(boson-boson) s-wave scattering length. In the following, 
we will always assume boson-boson repulsion, i.e. cibb > 
0. 

In the presence of a strong optical lattice and a suffi- 
ciently shallow external confinement in the x direction, 
so that at any lattice site its value can be considered con- 
stant, the bosonic field operators can be expanded in the 
basis of the single-particle Wannier wave-functions local- 
ized at each lattice site xi. Since the typical interaction 
energies involved are normally not strong enough in order 
to excite higher vibrational states, we can retain only the 
the lowest vibrational state in each lattice potential well 
(single-band approximation). In the case of stronger ex- 
ternal confinements, or interactions, one should include 
higher Bloch bands as well in the expansion of field op- 
erators, a case we do not consider in the present context. 
Moreover, as far as the harmonic trapping potential is 
concerned, we have shown in a previous work |27| that 
the introduction of higher energy levels of the harmonic 
oscillator does not modify the basic phenomenology of 
the system. Under these conditions, we can avoid work- 
ing directly with the exact Wannier wave functions and 
replace them, with an excellent degree of fidelity, with 
their harmonic-oscillator approximations at each optical 
lattice well. Then, the Wannier wave functions w(f) fac- 
torize in the product of harmonic oscillator states in each 
direction: 



= ^2diiv(x - Xi)w(y)w(z) , 



(7) 



moj 



(a-a + A 2 y 2 + A 2 z 2 ) 



(4) 



where Xi is the center of the i-th lattice well and di is 
the bosonic annihilation operator acting at the i-th lat- 
tice site. In each lattice potential well, the Wannier local 
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ground states are Gaussians in the harmonic approxima- 
tion: 



where 



U](x — Xi) 

w(y) 
w(z) 



1 



1 



1 



exp 



exp 



-{x - Xi) 2 



-y 



2 



exp 



2L 2 



2L 2 



(8) 



In Eq. JHJ we have introduced the harmonic oscilla- 
tor lengths of the groun d state in the y and z direc- 
tions Lj_ = y/h/ (m\uj) , and the oscillator length in 
the harmonic approximations of the periodic potential 
l x that reads, as a function of the lattice parameters, 

l x = (a 4 E R /(TT 4 V )) 1/4 , where E R = (nh) 2 /2a 2 m is the 
lattice recoil energy. In this paper we will consider the 
physical situation of very shallow trapping potentials, 
such that L x = Lj_vA ^> aM, with M denoting the 
total number of lattice sites. As a consequence, the local 
density approximation (LDA) can be applied. Therefore, 
when exploiting the expansion Eq. Q to map the full 
microscopic Hamiltonian Eq. (JTJ into its lattice version, 
we will discard all terms that are of order (aM/L x ) 2 or 
higher. Qualitatively, this means neglecting those nonlo- 
cal effects that are induced by the presence of the trap- 
ping potentials, such as site-dependent hopping terms. 
The latter can become important in regions of the lat- 
tice very far out of the central core of the harmonic trap. 
However, the typical experimental situations involve only 
that part of the lattice that lies well inside the central 
core of the slowly varying confining potential |^. We 
can then write down the translationally invariant lattice 
version of Hamiltonian Eq. in the form 



H 



(9) 



ho 



In Eq. © Ui t j t k,i is the two-body interaction strength 
that involves four sites of the lattice that depends on the 
relative distance between the sites involved. Recalling 
the expression of the two-body interaction Eq. @), to- 
gether with the form of the bosonic field operator Eq. (JJJ 
and of the lattice wave functions Eq. ©, one has 



U, 



U e 



7 /2 



(10) 



U = (2tt 



-S 9b b 



(11) 



is the local interaction strength, i.e. the amplitude of the 
interaction when i = j = k = /, 



e = exp {-a 2 /All 



is the lattice attenuation factor, and 



(12) 



7 = (i-j) 2 + (*-*0 2 + {i-l? 

+ (j - k) 2 + (j - I) 2 + (k- I) 2 , (13) 
is the "four-site distance" relative to all possible inde- 
pendent pairs of sites that can be chosen out of a set of 
four sites. It is worth noticing that e can be re-expressed 
in the form e = cxp(— 7t 2 y / s/4), i.e. in terms of an exper- 
imentally measurable and tunable quantity, the depth s 
of the lattice wells: s = Vo/Er. 

On the other hand, in Eq. Uj is the strength of 
the contributions of the kinetic and the external poten- 
tial terms to the energy. For i — j it gives rise to a 
constant zero point energy term that can be discarded 
by redefining the zero of the energy; if i ^ j it represents 
the probability amplitude for an atom to tunnel from 
the i-th lattice site to the j-th one along the x direction. 
Obviously, also this probability amplitude is a function 
of the distance between the involved sites, but it is im- 
possible to write for it a closed analytical formula like 
Eq. (|10fl . However, from the form of Eq. (JSJ it is easy to 
show that tij is still proportional to some positive power 
of the lattice attenuation factor e, with the exponent de- 
pending only on the distance between the sites. Hence, 
in general one can write 



Ju 



Ai-3) 2 



(14) 



where J\i-j\ decreases as a polynomial function of the 
modulus of the distance between the sites. Taking into 
account the form of the kinetic energy Eq. (J2J , the forms 
of the external potentials Eq. Eq. (@J, and Eq. (JSJ), 
together with the expression of the bosonic field operators 
Eq. Q and of the lattice wave functions Eq. JSJl, for 
\i — j\ — 1 (nearest neighbors) one has 



J 1 = V,(4-1- e-™<°') - ?!*!k - 2A. B - j"! c 2V - l) . (15) 



Writing the single particle hopping amplitude as in tions to the energy are function of the lattice attenuation 
Eq. ijTil) remarks the fact that the one-body contribu- factor e as well. In typical experimental situations, the 
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lattice spacing a is usually much larger than the local 
ground state length l x at each lattice site. Hence 
and this fact allows to exploit the lattice attenuation fac- 
tor as a meaningful dimcnsionlcss expansion parameter. 
The first two nontrivial contributions to the one-body 
part of the energy are proportional to e, corresponding 
t° \i ~ j\ — 1 an d e 4 , corresponding to \i — j\ — 2. The 
first contribution is the usual nearest-neighbor hopping 
of the standard Bose-Hubbard model, while the second 
one is a next-to-nearest-neighbor hopping term. 

Concerning the two-body contributions to the energy, 
the classification in powers of e looks in principle more 
complicated. However, it is easy to see that all the terms 



related to pairs of nearest-neighbor sites are proportional 
to powers e l at most of order 1 — 2. On the other hand, 
energy terms involving pairs of next-to-nearest-neighbors 
or triples of three adjacent sites are always smaller, be- 
cause the leading terms of this two classes of energy con- 
tributions are, respectively, proportional to e 6 and e 3 . 
Hence, in the lattice Hamiltonian description of interact- 
ing bosonic atoms in periodic optical potentials, we need 
to consider, in first approximation, only the energy terms 
proportional to e l with I < 2. At this order of approxi- 
mation, the lattice Hamiltonian Eq. reads, with the 
terms ordered in increasing powers of e, 



H = ^Y^hiihi - 1) - y£^(alai+i + H.c.) + U e§ ^ [(ofni&i+i +a|_ 1 rj. i a. i ) + H.c. 

i i 

+ 2U e 2 fkrk+i + U Q e 2 ^ (ijii+i + H., 



(16) 



where we have introduced both the on site occupation 
number operator hi — aja, and the on site pair annihi- 
lation operator A$ = &?. The first two terms of Eq. I|16|) 
represent the usual BH Hamiltonian with on site inter- 
action and single-atom nearest-neighbor hopping. The 
remaining terms give the corrections to this model, as- 
sociated to higher powers of the lattice attenuation fac- 
tor e. The term proportional to e 3 / 2 is the single-boson 
nearest-neighbor hopping conditioned by the on site oc- 
cupation; the first term proportional to e 2 is the density- 
density nearest-neighbor interaction; and, finally, the sec- 
ond term proportional to e 2 is the nearest-neighbor hop- 
ping of pairs of bosons. 



III. THE FREE ENERGY 

When dealing with systems of interacting bosons, it is 
convenient to work in the framework of the grand canon- 
ical ensemble 0, H 0, ^2^. Let us then introduce the 
grand canonical Hamiltonian K 

K = H-i±'Y j h l , (17) 

i 

where [i is the chemical potential needed to fix the aver- 
age total number of bosons in the lattice. All the summa- 
tions entering in Eq. (|17l) can be arranged in two differ- 
ent sets (Ki) and (Ki nt ). The first one, (Ki), contains all 
"local" terms that depend only on the on site occupation 
number operators hi and hi+±. The second one, (Kmt) 
contains all the "non local" hopping terms. According to 



this grouping, one can write 

k = k x + k nl . (is) 

We will first analyze the EBH model when the second, 
the third and the last terms of the right-hand side of 
Eq. (|16fl . i.e. the kinetic contributions to energy, may 
be neglected and only the local terms are retained (in a 
sense that will be clarified below). 

A. Local energy terms and mapping to a quantum 
Ising antiferromagnet 

Considering Eq. I|ltj|) and Eq. (|17|l . the local energy 
part in Eq. (jT%|l reads 

k = ^^hiiht - 1) - n hj + 2Uoe 2 ^2hih l+1 . 

i i i 

(19) 

Since in the "local" part of the Hamiltonian we include 
the nearest-neighbor interaction term, we should qualify 
that here, by "local" we mean all effects that do not 
involve particle exchange between sites. 

To fix techniques and notations, we first briefly recall 
how to determine the energy spectrum of the local part 
of the Hamiltonian in the standard BH model, i.e. when 
we neglect the term proportional to e 2 in Eq. (|19fl . In 
this case, it is well known that each term in the sum of on 
site interaction energies reaches its minimum for rii = n* 
with 

»--5 + &- (*) 
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This trivial observation naturally leads to introduce the 
complete orthonormal Fock basis and determines 

an obvious but important classification. If n* is close to 
an integer value, say no, we have that the energy gap 
A between the ground state |no) and the first excited 
state is of the order of the coupling constant Uq. This, 
incidentally, justifies neglecting, in first approximation, 
any correction proportional to any power / > 2 of e. In 
this situation, at zero temperature, the system is in a 
Mott-Insulator phase with exactly n atoms per site. On 
the other hand, if n* is closer to a half-integer value, 
then the energy gap A can be comparable with other 
contributions to the interaction energy, and one can write 



* 1 
11 = 2 



no 



2A 

Uq 



(21) 



where no is integer, |A|/[/o <C 1 and the two number 
states nearly degenerate in energy are |no) and |n + 1). 
This near degeneracy occurs in pairs: the states |no + 2) 
and |no — 1) (when it exists, i.e. when n > 1) are as 
well nearly degenerate and are separated from the pair 
{|n ), |n + 1)} by a gap of the order of Uo- The two 
nearly-degenerated ground states |no) and |no + 1) are 
separated from each other by an energetic distance equal 
to 2 |A|, while the gap between the two nearly-degenerate 
first excited states |no + 2) and |n — 1) is equal to 6 |A|. 
These results hold analogously for the pairs {|no + 3), 
|n — 2)} and {|n + 4), |n — 3)}, and so on, as long 
as the second element \uq — k) of each pair exists. From 
this classification, it emerges the fundamental role played 
by the two nearly degenerate states |no) and |no + 1). 
Moreover, the previous analysis allows to recast the local 
Hamiltonian Eq. I|19|) in a very useful form. Introducing 
the operator rhi = hi — (no + \) and fixing the zero of the 
energy, the local part of the grand canonical Hamiltonian 
reads 



o \ -» 



K miihi+i 



(22) 

where 6 = A - 2U e 2 (n + |), and K = 2U e 2 . 

The above discussion and Eq. (|22J) allow a clear un- 
derstanding of the local part of the EBH model at zero 



J 



temperature, showing that all the states different from 
the two quasi degenerate ground states |no) and |no + 1) 
do not contribute. Then, the first term in Eq. (|22|l can 
be neglected, and the remaining part of the local Hamil- 
tonian (|22|l describes an assembly of interacting two-level 
systems. Actually, because K > 0, it is mapped exactly 
in a spin 1/2 antiferromagnetic quantum Ising model in 
the presence of an external field —26: 



H eq = -26 J2 % + K^la, 



(23) 



At zero temperature, this model describes a system that 
undergoes a quantum phase transition at the critical field 
6 C = K/2. The ferromagnetic phase \6\ > 6 C corre- 
sponds to the Mott-Insulator phase with the same, con- 
stant number of atoms no on each lattice site. The an- 
tiferromagnetic phase 1 6 \ < 5 C corresponds to a density 
wave phase with no atoms on a site and no + 1 on its 
neighbor. Analogously to the spin system in the anti- 
ferromagnetic phase, the optical lattice for the bosonic 
atoms in the density wave phase is divided in two sub- 
lattices of "staggered" atomic densities n and n + 1 . In 
the following, we will denote the two phases, respectively 
by PMI (Pure Mott-Insulator) and by DWMI (Density 
Wave Mott-Insulator). 



B. Nonlocal energy terms, ferromagnetic- and 
antiferromagnetic-like models, and the mean field 
free energy 

When the nonlocal hopping terms are reintroduced, 
the tensor product states of the local occupation number 
states (local Fock states) are no longer eigenstates of the 
total Hamiltonian. The true eigenstates cannot be deter- 
mined analytically, and consistent approximations must 
be envisaged to approach the problem in the product ba- 
sis of the local states. Let us first rewrite the nonlocal 
terms appearing in the grand canonical Hamiltonian 1)18(1 
in terms of the on site "magnetization" operators rhf. 



fai+i + H.c)j + U Q ei («!(™* + rh i+1 )a l+1 + ff.c.) + y ^ ( A l A *+i + Hx 

i i i 

i 

where the "dressed" hopping amplitude J reads From Eq. (|25Jl we see that the density-dependent part 

of the hopping amplitude gives a negative contribution 
if the boson-boson interactions are repulsive (Uq > 0). 



it. 



(24) 



J = e 



Ji - 2U a e~- ( n + - 



(25) 
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Stability of the ground state energy thus requires 

Ji > 2U ei (n + . (26) 

Typically, such a stability requirement can be easily sat- 
isfied in most experimental situations, unless one goes to 
very large occupation numbers no and very strong inter- 
action couplings Uq. Hence, in the following discussions 
and examples, we will always consider situations in which 
the stability condition Eq. (|26|) is satisfied. 

Before introducing mean field approximations, we first 
need to deal with the second term in Eq. (|24[) . the condi- 
tioned hopping term. This can be done, in a Bogoliubov- 
like framework, by replacing the operator TOj with its 
average value Xi = < t^H >; thus neglecting quantum 
fluctuations. This is justified in the situation in which 
the magnitude of the on-site interaction amplitude Uq is 
sufficiently large that the probability of finding on site 
occupation numbers that do not fall in the ranges iden- 
tified by the pairs {no, no + 1} and {no + 2, no — 1} is 
negligible. This is the physical situation that one usually 
meets in realistic experimental conditions. In implement- 
ing this approximation we must thus distinguish between 



two different instances, according to the previous discus- 
sion on the local terms. 

/-If one has \8\ > K/2, then the approximate model 
describing the system is ferromagnetic-like and, in the 
limit of vanishing nonlocal hopping terms, the associated 
ground state reduces exactly to the PMI phase. In the 
ferromagnetic-like model, the expectation value of the oc- 
cupation number Xi is constant on all sites of the optical 
lattice: Xi = X 

II- If one has \S\ < K/2, then the approximate model 
describing the system is antiferromagnetic-like and, in 
the limit of vanishing nonlocal hopping terms, the as- 
sociated ground state reduces exactly to the DWMI 
phase. In the antiferromagnetic-like model, the expec- 
tation value of the on site occupation number takes op- 
posite values on neighboring sites of the optical lattice: 
Xi = X-, Xi+i = ~~ X an d henceforth the conditioned hop- 
ping term in Eq. I|24|) vanishes. 

Then, in compact notation, the grand canonical Hamil- 
tonian K F for the ferromagnetic-like case and the grand 
canonical Hamiltonian K A for the antiferromagnetic-like 
case read: 



K 



F.A 



E 



T F,A 



H+l 



£ (otoi+i + H.c) + | (A 1 A+i + H. 



(27) 



Here, J F,A are the single particle hopping amplitude of 
the two models. When \8\ > S c , one must take the choice 
J F = J — AUqe^x- When \8\ < 8 C , one must take the 
choice J A = J. 

Before proceeding, we would like to make the following 
side observation. The grand canonical operator Eq. I|27|) 
takes into account both the local and the nonlocal parts 
of energy. We have seen that when we can neglect the 
nonlocal terms (i.e. in the absence of the kinetic terms), 
the local part of the Hamiltonian Eq. (|19|l is mapped in a 
spin-1/2 quantum Ising model Eq. We can exploit 

a similar mapping also for the total Hamiltonian Eq. I|27|) 
in the particular limit when the on site interactions are 
strong enough that the only states that contribute are 
those with no and no + 1 bosons per lattice site. In this 
special limit, the nearest neighbor atomic pair hopping 
operator A[Ai + i does not produce any effect. Hence, by 
the same identifications d\ — &\ = of + ib\ and mj = of 
that map Hamiltonian Ki Eq. i|19|) to the quantum Ising 
model Eq. J23J), the two Hamiltonians K F > A Eq. J23 are 
mapped in the two spin-i XX Z Hamiltonians 

Hx'xz = -26Y,of + Kj2 *t°! + i 

i i 



-J F ' A X)(fiffifn + W +1 ). (28) 

i 

The above limiting mapping of Bose-Hubbard models in 
XXZ Hamiltonians in external field has been investi- 
gated extensively by van Otterlo et al. |49|. who pre- 
dicted the existence of a supersolid phase in two dimen- 
sions. 

We are now in a position that allows to introduce 
mean field approximations on the generic terms contain- 
ing pairs of operators an adjacent sites, namely: a|Sj+i; 
A|j4j+i; and miihi+i. For the latter term, we must make 
a bookkeeping for the two different model Hamiltonians 
corresponding to \8\ > \8 C \ and \8\ < \8 C \. In the first case 
we must consider x =< ™; > Vi. In the second case 
the order parameter has opposite signs on adjacent sites. 
For the first two pairs of terms, in order to implement 
correctly the mean field approximation, we must make 
sure that concavity of the energy holds, guaranteeing that 
the extremal conditions correspond to a true minimum of 
the energy and not to a maximum. For the single parti- 
cle hopping, recalling that —J Ft < we must consider a 
uniform order parameter < a. L >=< a I >= <f> Vi. For the 
pair hopping term we have that its amplitude K is always 
positive; hence, in order to obtain the right concavity, we 
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need to choose an order parameter that takes opposite 
signs on adjacent sites: < A, t >= tp, < A i+ i >= —ijj. 
We have chosen to restrict to real order parameters even 
if, due to the non-hcrmiticity of the involved operators, 
in principle complex order parameters would be allowed. 
Obviously, the imaginary parts of the order parameters 
may be easily taken into account. However, we have ver- 
ified that even in these cases the extremal conditions al- 
ways lead to real results. For this reason we can restrict 



our analysis right from the start to real order parameters, 
a situation that is in complete analogy with the one en- 
countered in the study of the standard on site BH model 

sua 

We can now write down the mean field expressions for 
the ferromagnetic-like and antiferromagnetic-like grand 
canonical total Hamiltonians (with M being the total 
number of lattice sites) : 



K 



F,A 



2 

K 



E N - j) - jF ' A E ( & t + **) * - v - **) E ^ - (* * k X ) e 



m, 



iesi 



ieS2 



M 



(29) 



iGSl 



it.S'2 



where the minus sign holds for the ferromagnetic and 
the plus sign for the antiferromagnetic case. By SI and 
S2 we denote the two different sublattices in which the 
original lattice is split with regard to the if) and x order 
parameters (for the latter, only in the case 2\S\ < K). 

The grand canonical total Hamiltonians are written 
down as sums of local on site energy terms, and the order 
parameters {%, <f>, ip} must be evaluated self-consistcntly. 
In principle, the Fock spaces associated to the on site 
occupation numbers are infinite-dimensional. However, 
the leading term in Eq. (|29|l is the one proportional to 
Uq , so that all number states with eigenvalue greater than 
Uq can be neglected, leading to consider only the set of 



the four lowest lying states that include the two nearly 
degenerate local states |n > and \n + 1 >, and the two 
nearly degenerate local states \no + 2 > and |n — 1 > 
(when the latter exists, i.e. when no > 1). 

Starting from the two grand canonical Hamiltonians 
Eq. (|29|l we can evaluate analytically the two free ener- 
gies F of the system at any inverse of the temperature 
(3 = (ksT)^ 1 cither in the ferromagnetic-like or in the 
antiferromagnetic-like case. One elegant technique to do 
so is to resort to the resolvent approach, as illustrated in 
the Appendix. Considering the free energy per site f F ' 
in the thermodynamic limit, one has: 



(J F ^) 2 (n + l) + gy(n + l) s 
U 



1 

2^5> g 

r— 1 



2 cosh /3(A r + -f) 



(30) 



where 



A, = ^(5-K X r) 2 + (J F < A 0) 2 (no + l) 



= - — { 2(J F ^) 2 KMno + l) 2 + [(J F ' A <j>) 2 - K 2 tf(no + 1)](S - K Xr ) 



(31) 



In Eq. (|30|l the index r runs over the two sub-lattices 
Si and S2 in which the original lattice is split. In the 
ferromagnetic case 2\6\ > K the two sublattices coincide 



and Ai = A2 = A, ol\ = oli = a, where 



(32) 
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a = -j{2(J F <j)) 2 K^(n + l) 2 + 

((J F 4) 2 ~K 2 ^(n Q + l))(5~K X )} . (33) 

The free energy per site so obtained depends, obviously, 
on the three order parameters </>, ip and \, that must be 
evaluated self-consistently. Regarding <fi and ip, this is an 
easy task; it is accomplished by simply determining the 
minimum of the free energy in each case. The existence 
of the minimum is assured by the right concavity of the 
free energy and hence it is enough to impose Of /d<j) = 
and df/dip = 0, in order to determine their extremal 
values. On the contrary, the order parameter \, that de- 
pends both on </> and ip: x — x(4 ) ->' l l > )> cannot be simply 
evaluated by fixing the extremality conditions. One must 
instead resort to its definition, and solve analytically for 
it, i.e., we must use the fact that x is defined as the aver- 
age value of m, and exploit this definition to determine it. 
This evaluation may be performed using the same mathe- 
matical technique employed for the evaluation of the free 
energy (see the Appendix). The self-consistent equation 
so obtained together with the extremal condition with re- 
spect to the hopping order parameter are the set of rela- 
tions that are needed to analyze the phase diagram of the 
system. As we have already mentioned, the expressions 
derived for the Hamiltonians and the free energies are 
obtained and are valid in the moderately strong coupling 
regime, where only the first four lowest lying states are 
considered (see the Appendix for more details) . This con- 
dition is consistently met when the ratio w of the dressed 
hopping to the on site interaction coupling strength does 
not exceed unity: w = J(no + l)/C^o < 1- 



IV. RESULTS AND DISCUSSION 

The study of the different possible solutions of the 
three equation needed to determine the different order 
parameters supplies the information needed about the 
phase diagram of the system. Obviously, it is not pos- 
sible to follow analytically all the solutions as functions 
both of the temperature and of the Hamiltonian param- 
eters, and exact numerical solutions will be used to track 
the phase diagram in the whole range of physical param- 
eters. 



A. Phase diagram: qualitative aspects, absence of 
pair superfluidity, and the role of many-body 
interactions 

As it is well known, the standard BH model sustains a 
phase transition between a single-boson superfluid phase 
and a normal (disordered) phase that at vanis hing tem- 
perature reduces to a Mott insulator phase H,0j0|. The 
first issue we wish to address here is whether, due to the 
presence of energy terms corresponding to the hopping 
of pairs of atoms between adjacent sites, the EBH model 



can sustain a new superfluid phase characterized by a non 
vanishing value of the pair-superfluidity order parameter 
ip either in the absence (<fi = 0) and/or in the presence 
(4> =/= 0) of the standard superfluidity of individual atoms. 
In fact, we find that within the EBH setting this is never 
the case, and it is always ip = 0. This negative result 
can be easily understood by looking at the condition of 
extremality obtained by differentiating Eq. I|3U|I with re- 
spect to the pair-superfluidity order parameter for each 
sub-lattice. By recalling that within the two sub-lattices 
(r = 1, 2), tpx = —tp2 = ip, one finds: 

, K(n + l) 2 , 1 da r , / a r . \ 

U AKUo ^-J dip V u o J 

(34) 

We immediately observe that the left-hand side of 
Eq. I|34l) involves both terms proportional to the zero- 
order power of Uo and to the inverse of Uo, while in 
the right-hand side only the term proportional to the 
inverse of Uo appears. Since our analysis is carried out 
in the strong-coupling limit (large Uo) and the hyper- 
bolic tangent takes values in the range [—1,1], the ex- 
tremality condition Eq. (|34|l is effectively of the form 
(1 — A)ip = Bip. This relation, regarded as a linear alge- 
braic equation for the pair-superfluidity order parameter 
ip, is satisfied only if ip vanishes, due to the fact that the 
real coefficients A and B are both very small: A, B <C 1 
(In fact, in most situations it is even 4,6 « 0.1). To 
illustrate how this takes place, we focus on the case in 
which the single particle order parameter <p vanishes (it is 
easy to check numerically that the same conclusions hold 
true for </> 7^ as well). For null single particle superflu- 
idity, (p = 0, by keeping in mind that K — 2Uoe 2 , one has 
A = 2e 2 (n + l) 2 and B = e 2 (n + l) tanh(/3(A+ ^))U=0- 
Let us fix, for instance, no = 9, and evaluate the co- 
efficients A and B at two different values of the lat- 
tice attenuation parameter, e = 0.01 and e = 0.001. 
Then, in the first case we have A = 0.02 and B = 
0.001 tanh(/3(A + ^))|^ = o. In the second case we have 
A = 0.0002 and B = 0.00001 tanh(/3(A + #-))U=o- 

The circumstance according to which the pair- 
superfluidity order parameter vanishes holds in general. 
In fact, one can show that taking into account corrections 
proportional to any power I > 2 of e, leads to new hop- 
ping terms that are different from the single-particle and 
pair hoppings that we have considered so far. An exam- 
ple of such hopping terms of higher order is provided by 
the operator describing the collective tunneling from, say, 
site i to site j of pairs consisting of two atoms localized on 
nearest neighbor sites, and so on. Each of these new hop- 
ping terms is associated to a suitable order parameter to 
be determined self-consistently. For these order parame- 
ters, the same arguments exploited for ip hold, implying 
that the extremality conditions are always and only sat- 
isfied if all the order parameters for the superfluidity of 
composite particles vanish identically. This result leads 
to a behavior, with respect to superfluidity, that is ruled 
by the single-atom superfluid order parameter 0, and is 
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thus qualitatively similar to that exhibited by the stan- 
dard BH model at finite temperature, under the effect of 
a superimposed harmonic confinement |27j . 

From a fundamental physical point of view, the impos- 
sibility to obtain a superfluid phase in which the tunnel- 
ing particles are composite bosons made up by two or 
more bosonic atoms stems from the fact that the local 
eigenstates of lowest on-site energy are always the two 
consecutive Fock states |no) and \no + 1). This always 
makes the hopping of any aggregate of atoms energeti- 
cally unfavorable. To engineer a superfluid phase of com- 
posite bosons is thus necessary to overcome this limita- 
tion and realize a situation in which the two lowest local 
eigenstates of lowest on site energy are \no) and \no + fc), 
with k > 1 being the dimensionality of the generic com- 
posite. To achieve this goal is then necessary to engineer 
and take into account many-body interactions compara- 
ble in magnitude to the standard bilinear ones (two-body 
collisions) that are usually the only interactions included 
in the description of dilute systems of interacting bosons. 

In the following, we will first consider the ferromag- 
netic model to track the SF-PMI quantum phase tran- 
sition as the zero-temperature limit of the the finite- 
temperature transition between the disordered and the 



J 



superfluid phase. Later on, we will consider the anti- 
ferromagnetic model to determine the SF-DWMI phase 
transition, and, finally, we will analyze the full quantum 
phase diagram of the EBH model at zero temperature. 



B. Finite- and zero-temperature transitions to 
superfluidity 



Starting from the "ferromagnetic" grand canonical free 
energy, the critical diagram for single atom superfluidity 
is reported in Fig. ^ for different sets of values of the 
Hamiltonian parameters. From Fig.^ we observe a low- 
ering of the critical temperature T C F with increasing am- 
plitude of the energy gap A. From a physical point view, 
this situation is due to the fact that for high values of A, 
bosons experience a high potential barrier that contrasts 
the hopping from a site to its nearest neighbor with a 
consequently increasing difficulty for the whole system 
to go toward an ordered phase and hence the superfluid 
transition occurs in "colder" zones. Solving the equation 
df /d<j)\^ = o = with x evaluated at ij) = <fr = 0, we get 



1 



8-H X (P,0) 



tanh 



2/(J-ifx(0,0))-2(n + l)/C/o 



{n + l)/(S-K X (0,0))-2/U 



(35) 



r 



where H = 2Uo£% ■ In Fig. [2] we show the behavior of the 
critical temperature in units of the lattice recoil energy 
Er as a function of the filling factor n = no + 1/2 + x 
for different values of the Hamiltonian parameters. Due 
to the existence of a region with density wave order, the 
condition x = 0orrt = n + l/2is verified throughout 
an entire (although small) region in the space of param- 
eters rather than at a given point in it. In this region 
the critical temperature, as we will see in the following, 
becomes function of the DWMI order parameter x, while 
the filling factor remains constant. Hence Fig. |3 repro- 
duces the behavior of the critical temperature only in the 
ferromagnetic-like instance, and the value of the critical 
temperature with semi-integer value of the filling factor 
requires a longer analysis that will be presented in the 
following subsection. Fig. [21 shows the competition be- 
tween thermal effects and ordered mobility. At fixed n, a 
larger hopping amplitude corresponds to a higher critical 
temperature. In Fig.[3|we show instead the behavior of 
the critical temperature as a function of the optical lat- 
tice depth s = Vo /Eji for different values of the chemical 
potential. As the depth of the lattice increases, hopping 
and mobility are suppressed, and the critical temperature 
of the superfluid transition lowers. 



Looking back at Fig. we must notice that, obviously, 
the critical temperature vanishes for integer value of the 
filling factor (no superfluidity allowed). Requiring in- 
stead that assumes an infinite value in Eq. (|33|l . we 
obtain the critical condition on the local gap or "exter- 
nal magnetic field" 8 for the zero-tempcrature quantum 
phase transition from a PMI phase to superfluidity: 



5* = 



(J-H/2)(n + l) 



2 1- 



ng(J-H/2) 
Uo 



K 
~2 



(36) 



If the local "gap" 5 is smaller than 8f , the system is in a 
superfluid phase; otherwise, a Mott insulator is realized. 
Clearly, this result holds provided that S > i.e the 
system cannot access the density wave region. The be- 
havior of 6f as function of the depth of the optical lattice 
is showed in Fig. for different values of the anisotropy 
parameter occurring in the external harmonic confine- 
ment. Concerning Fig. |3| we observe that the functional 
behavior is not plotted down to small values of the lattice 
depth parameter s. This is due to the fact that for values 
of s in the approximate range [0, 15], the weak coupling 
ratio w = J(jiq + 1)/Uq may exceed unity, so that the 
strong coupling approximation breaks down. The graph 
has thus been plotted in the interval of values of s such 
that < w < 0.6. We observe that the critical temper- 
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FIG. 1: Single atom superfluid order parameter <j>, as a func- 
tion of the dimensionless "ferromagnetic" critical temperature 
ksT F /En rescaled in units of the lattice recoil energy En, for 
J = 0.01f/o- From top to bottom, behavior for A = (black 
solid line), A = 0.009(7 (dashed line), A = 0.02[/ (dotted 
line), and A = 0.04(7o (solid line). Notice, as expected, that 
the critical temperature lowers as the energy gap A increases. 



FIG. 3: Dimensionless critical temperature ksT^ /En, 
rescaled in units of the lattice recoil energy En, for the tran- 
sition to superfluidity as a function of the dimensionless lat- 
tice depth s = Vo/En- From top to bottom, behavior for 
A = (solid black line), A = 0.0005(7 (solid gray line), and 
A = 0.0009C/ (dashed black line). 
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FIG. 2: Dimensionless critical temperature IcbT^ /En, 
rescaled in units of the lattice recoil energy En, for the tran- 
sition between the disordered (high-temperature) phase and 
the superfluid (low-temperature) phase as a function of the 
filling factor n. From top to bottom, functional behavior 
for J = O.Olt/o (dashed line), J = 0.007t/ (dotted line), 
J = 0.004£/ (dashed-dotted line), and J = 0.001[/o (solid 
line). As the overall hopping J increases, the critical temper- 
ature rises. 



FIG. 4: Dimensionless critical field 8^ / En, as a function of 
the lattice depth a. From left to right, behavior for A = 65 
(dashed black line), A = 39 (solid gray line), A = 13 (solid 
black line). 

< w < 0.6 holds. We see that the range of permissible 
values of s grows with increasing anisotropy. 



ature decreases for increasing s. This is due to the fact 
that the greater the lattice depth, the more the on-site 
interaction tends to dominate on the hopping. Hence, 
in order to achieve the onset of the transition to the or- 
dered superfluid phase it is necessary to operate at lower 
temperatures. Moreover, once the lattice depth is fixed, 
the critical temperature is lower for larger energy gap A, 
in agreement with the results presented in Fig. ^ Con- 
cerning Fig. 0] each of the three curves, corresponding 
to a different value of the transverse trapping frequency 
(anisotropy parameter), is plotted for a different range 
of the lattice depth s. Each of these different ranges 
corresponds to to the different zones in which, for the 
various values of the anisotropy parameter, the relation 



C. Unified finite-temperature phase diagram 

The transition from an ordered superfluid phase to a 
Density Wave Mott Insulator can be determined start- 
ing from the antiferromagnetic grand canonical free en- 
ergy along the same lines followed to analyze the SF- 
PMI phase transition in the ferromagnetic grand canoni- 
cal setting. Hence, we will analyze the zero-temperature 
SF-DWMI quantum phase transition by first determin- 
ing the "antiferromagnetic" critical temperature and 
critical field, 5^, the associated finite-temperature phase 
diagram, and by finally taking the zero-temperature 
limit. Concerning the first step, straightforward evalu- 
ation yields: 
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1 - 



J(n + 1) 
U 



= J ^ I tanh /3 C ((5 + K Xr (0, 0))] , ) = 



(37) 



for the inverse (3^ of the "antiferromagnetic" critical 
temperature, and 



J 



JU (n + 1) + y/4K*(J (1 + n ) - £/ ) 2 + (J E/o(rc + 1))- 
4(U - J (1 + n )) 



(38) 



r 



for the "antiferromagnetic" critical field, that differs an- 
alytically from the "ferromagnetic" one expressed by 
Eq. H3(jl) . Of course, in those regions of the space of pa- 
rameters that allow to neglect the nonlocal energy terms, 
5£ and coincide exactly. The relations expressed by 
Eq. i|37|) and Eq. together with the corresponding 
ones previously obtained for the SF-PMI phase transi- 
tion allow us to construct the full phase diagram of the 
EBH model both at finite and at zero temperature. 

Concerning the finite-temperature scenario, we report 
in Fig. [3 the behavior of the critical temperature for the 
transition from the high-temperature disordered phase to 
the low-temperature ordered superfluid phase as a func- 
tion of the energy gap A for different values of the lattice 
parameters. Recalling the relation S = A — 2U~oe 2 (n + ~) 
that connects the local external field with the energy 
gap, we can follow the entire evolution of the critical 
temperature as a function of the external field, mov- 
ing smoothly through the ferromagnetic and the anti- 
ferromagnetic regimes. The thermodynamic evolution of 
the system may be considered made up of three con- 
tinuous intervals. Going from left to right on the ab- 
scissa in Fig- El the interval of negative values of A, that 
maps in the interval of negative values 6 < —K/2, corre- 
sponds to a ferromagnetic critical temperature T c = Tf . 
The central part of the interval, around A = 0, corre- 
sponds to the interval of negative and positive values 
—K/2 < 5 < K/2 and to an antiferromagnetic critical 
temperature T c = T^. Finally, the right part of the in- 
terval of positive values of A that maps in the interval 
of positive values 5 > K/2, corresponds again to a ferro- 
magnetic critical temperature T c = T C F . 

In the central interval, the finite-temperature analogue 
of the zero-temperature SF-DWMI quantum phase tran- 
sition is realized (antiferromagnetic- like coupling). In the 
two external regions the finite-temperature analogue of 
the zero-temperature SF-PMI quantum phase transition 
is realized (ferromagnetic-like coupling). In the former, 
antiferromagnetic-like case, the behavior of the critical 



temperature as the gap varies in the range [— K/2, K/2] 
should be represented by a flat, constant line in the cen- 
tral region of Fig. \5\ joining the two curves represent- 
ing the critical temperature in the two ferromagnetic- like 
external regions. However, due to its extremely small 
extension, this "antiferromagnetic connection" appears 
shrunk to a single point, the overall maximum of the 
critical temperature. Hence, the only visible landscape 
in the regions above the critical temperature in Fig. [S] 
is the one relative to the finite-temperature analogue of 
the zero-temperature arrangements in which a Pure Mott 
Insulator arrangement is favored. In these regions, the 
critical temperature lowers as the modulus of the gap in- 
creases. This fact is in agreement with the considerations 
developed for Fig.^ In particular, when the modulus of 
the gap is larger than the value of the critical field for the 
SF-PMI quantum phase transition, the critical tempera- 
ture vanishes and the lattice is characterized by the same, 
constant integer filling factor on all sites. Moreover, by 
using the same arguments presented in the previous sub- 
section for the behavior of the critical temperature as 
a function of the filling factor, we can understand the 
lowering of the critical temperature as a consequence of 
the lowering of the hopping amplitude from a site to its 
nearest neighbor. 



D. Unified zero-temperature quantum phase 
diagram: Superfluid, Pure Mott Insulator, and 
Density Wave Mott Insulator phases 

To conclude our study, we can now consider the full 
diagram of quantum phases at zero temperature, by tak- 
ing the limit (3 C — > 0. At zero temperature, the quantities 
that determine the transition from a kind of ordering to 
another one are the Hamiltonian parameters, that are 
controllable quantities. When the ratios of these control 
parameters are suitably tuned, macroscopic changes take 
place in the ground state of the system. These changes 
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FIG. 5: The critical temperature ksTc / ' Er for the transition 
from a disordered phase to the ordered superfluid phase as 
a function of the local energy gap A/Er. From top to bot- 
tom, behavior for J = 0.01[/o (black solid line), J — 0.009£/o 
(dotted line), and J = 0.007?7o (gray solid line). 



give rise to the zero temperature phase diagram that we 
report in Fig. H3 Here the control parameters are the 
magnitude of the nearest neighbor hopping amplitude J 
and the local energy gap 8. In establishing the boundaries 
between the different phases, we must take into account 
that 8% and 5£ essentially coincide in a wide range of 
values of the lattice parameters. These two quantities 
determine, respectively, the boundary lines at the quan- 
tum phase transitions from the SF to the PMI phase, 
and from the SF to the DWMI phase. The quantity 8 C 
instead determines the boundary line at the quantum 
phase transition from the DWMI to the PMI ordering. 
An important novelty emerges with respect to the phase 
diagram of the standard BH model. In fact, in this last 
case there exists only one boundary line, the one sep- 
arating the SF from the PMI phase. However, in the 
quantum phase diagram of the EBH model, two further 
boundary lines appear. 

The first one is the coexistence curve for the SF and 
the DWMI orderings; the second one is the coexistence 
curve for the two insulating phases, the PMI and the 
DWMI. The three different boundary lines cross at two 
triple points where all the three phases coexist. From 
Fig. EJand Fig.0 describing the zero-temperature phase 
diagram for two different values of the lattice attenuation 
factor e, we of course see that the zone in which the 
system is in a DWMI phase is extremely small compared 
to the regions occupied by the SF and PMI phases. This 
could be already expected from the shrinking to a point 
of the corresponding antiferromagnetic-like plateau in the 
finite-temperature diagram reported in Fig. 

Comparison of the two phase diagrams reported in 
Fig. and Fig. [7| shows that the lattice attenuation fac- 
tor plays a crucial role concerning the area in the space 
of parameters in which the system is in a DWMI phase. 
The smaller is the value of e the smaller and less observ- 
able becomes the region in which the DWMI phase takes 
place. This fact implies that the experimental observa- 
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FIG. 6: The zero-temperature quantum phase diagram of the 
EBH model in the strong coupling regime. Horizontal axis: 
dimensionless gap S/Uo. Vertical axis: dimensionless normal- 
ized hopping amplitude J/Uq. The vertical lines are the sepa- 
ration lines between the DWMI and PMI phases. The oblique 
lines are the separation lines between the SF and PMI phases 
and between the SF and DWMI phases. Symmetrically placed 
on the sides of the cusp point are the two tricritical points. 
The phase diagram is plotted for a value of the lattice atten- 
uation factor e = 0.07 and on site occupation no = 9. Here 
the label PM stands for Pure Mott, DW for Density Wave, 
and SF for SuperFluid. 
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FIG. 7: The zero-temperature quantum phase diagram of the 
EBH model in the strong-coupling regime, plotted for e = 0.01 
and no = 9. Labels denoting the various quantum phases have 
been omitted, as the meaning of the oblique and vertical lines 
is the same as in Fig. |S] Notice, in particular, the dramatic 
shrinking of the density wave phase for a lower value of e, 
compared to the one fixed in Fig. 



tion of such a phase will require significant advances in 
the manipulation and control of systems of interacting 
bosons in optical lattice potentials. In particular, it will 
be important to combine optical lattices potentials and 
magnetic Feshbach resonances to enhance the on site in- 
teractions to strong coupling limits, while at the same 
time keeping the lattice attenuation factor in a range of 
not too exceedingly small values. 
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V. CONCLUSIONS AND OUTLOOK 

We have studied systems of ultracold spin-zero neu- 
tral bosonic atoms with repulsive interactions, harmoni- 
cally trapped and regularly arranged by means of a peri- 
odic optical lattice potential. Taking into account the se- 
ries expansion of the amplitudes of the interaction terms 
in powers of the lattice potential parameters and of the 
lattice attenuation factor, we have mapped the second- 
quantized total Hamiltonian in a new, specific form of 
Extended Bose Hubbard (EBH) Hamiltonian. We have 
then established various mappings of this atomic EBH 
model to models of interacting spin-^ systems. By using 
such a correspondence, we have analyzed in a unified way 
the Density Wave and Pure Mott Insulator phases sup- 
ported by the model, in analogy with the unified mean 
field treatment of ferro- and anti-ferromagnetism. 

We have developed the mean field theory description 
of the EBH model both at finite and zero temperature, 
determining the free energy density, and analyzing the 
finite-temperature behavior of the model, determining 
the phase boundaries between the ordered superfluid and 
the disordered high-temperature phase. We have demon- 
strated the theoretical possibility for two different tran- 
sitions to superfluidity within the EBH model, one due 
to the hopping of single atoms, and the other due to the 
hopping of atomic pairs. In fact, we have given a thermo- 
dynamical proof that only the first mechanism is realized 
if one truncates the expansion in the lattice attenuation 
parameter at lowest order. Finally, we have determined 
the zero temperature phase diagram of the EBH model, 
showing the existence of a new quantum phase, the Den- 
sity Wave Mott Insulator, which is not allowed within the 
framework of the standard BH model, and we have deter- 
mined the range of lattice and Hamiltonian parameters 
for which such a phase can be detected. The two differ- 
ent forms of localized phases, Pure Mott Insulator and 
Density Wave Mott Insulator, manifest themselves, re- 
spectively, in the different behavior of the atomic density 
in the lattice. The PMI phase is characterized, as usual, 
by the same, constant integer filling factor throughout 
the entire lattice; the DWMI is instead characterized by 
two different integer filling factors in two sublattices, say 
no for half of the lattice sites, and no + 1 on their neigh- 
bors (checkerboard phase) . We have studied the behavior 
of typical physical quantities of the system, illustrating 
how the different control parameters involved compete in 
determining the evolution of the system. 

Regarding future perspectives, it is to be expected that 
by taking the expansion of the second-quantized total 
Hamiltonian further up to higher powers in the lattice 
attenuation factor, a new, and accordingly very small re- 
gion in the phase diagram will emerge where pair super- 
fluidity, absent both in the BH and in the EBH model, 
can occur, at least at extremely high filling factors, as 
well as new types of intermediate range interactions and 
tunneling mechanisms. The same framework introduced 
in this paper may be extended to the case of systems 



of interacting bosons when the excited harmonic levels 
of the trapping potential |53 and/or the hi gher -excited 
Bloch bands of the optical lattice potential |45l |4?| are 
taken into account. It would be an interesting challenge 
to further extend the scheme developed in the present 
work for pure single-flavor bosons with repulsive interac- 
tions to the case of multi- flavor bosons and/or attractive 
interactions; to mixtures of bosonic and fermionic atoms 
interacting on a lattice 

El El EJ E2; and, finally, to 
the case of disordered and/or random optical lattices that 
allow for the study of disordered ultracold atomic gases 

mm- 

APPENDIX 

In thi s ap pendix, we will briefly review the resolvent 
method |55| needed to obtain the expression Eq. (|3L)[I for 
the free energy of the EBH model. 

In order to study the thermodynamic properties of the 
system in the grand canonical ensemble, we must deter- 
mine the corresponding partition function , Z: 

Z = Tr[exp{-f3K F ' A )} , (A.l) 

where K F ' A is given by Eq. (gU and (3 = l/k B T with 
ks the Boltzmann constant. By writing down the ex- 
plicit form of K F,A , the grand canonical partition func- 
tion reads 

Z = Tr{exp[- (3^2(11, + J F ' A cl) 2 

i 

+K^ T K X 2 )}}. (A.2) 

In the last equation hi represents the action of the oper- 
ator 

h = Y,k =Y,ChL + hj) (A.3) 

i i 

on the i-th lattice site. As we may deduce from Eq. I|29p. 
the first operator appearing inside the sum in the right- 
hand side of Eq. (|A.3|) is the Hamiltonian whose eigen- 
states are tensor products of local Fock states \uq + k > 
with k integer or zero: 

zZ hL = -£ (♦"? - 1 ) - 6 ( ™* + J2 ™> ) ' 

i i ^ ' \ieSl ieS2 / 

(A.4) 

while the second operator is the mean-field "decoupled 
version" of operators representative of the single-boson 
hopping, atomic-pair hopping, and of the density-density 
interaction rhiihj, respectively: 

i i 

+ K X Y™i- (TKx) ^2 rhi 
iesi ieS2 
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A" 



(A.5) 



In Eq. (|A.4|> and Eq. I|A.5|) . the indexes 51 and 52 de- 
notes the two sub-lattices in which the whole lattice is 
split. The meaning of the parameters J F,A and K have 
been already explained in Sections II and III. 
From the grand canonical partition function, the expres- 
sion for the free energy F F,A — —jjlnZ of the system 
is readily deduced. This thermodynamic potential will 
depend, in general, on the mean field parameters of the 
theory. 

We describe our system in the complete basis of the num- 
ber states. Keeping in mind that we are analyzing the 
physics of our system in the mean field approximation, 
correlations between different lattice sites are neglected 
and hence, the partition function Z of the systems fac- 
torizes into the product of M independent partition func- 
tions, each of these evaluated for a single site. In the Fock 
states basis and in mean field approximation framework, 
the grand canonical partition function then reads 



< n\ exp ( - 0{hi + J F ' A cj) 2 + Ktp 



TK X 2 ))\n> 



M 



(A.6) 



where M is the total number lattice sites, the sum is 
in principle performed over all Fock states, and it is in- 
tended that the thermodynamic limit must be eventu- 
ally taken. However, as already discussed in Section III, 
for sufficiently strong coupling we may limit ourselves 
to consider the four Fock states of lowest energy |no >, 
| no + 1 >, | no — 1 >, and \no + 2 > (actually, in the ultra- 
strong coupling re gim e, it is enough to consider only the 
two lowest states H^)- This choice is fully justified as 
long as the weak coupling parameter w = J(uq + 1)/Uq 
does not approach or exceed unity. In this way we can 
determine the expressions of the physical quantities of in- 
terest in the EBH model at first order in powers of 1/Uq- 
To evaluate the free energy, one needs to compute the 
trace of the operator exp(— (3hi). However, rather than 
diagonalizing hi in the space spanned by the four lowest- 
lying Fock states, it is more convenient to write the free 
energy per site f F,A in the following way: 



J 



cF,A 



M 



" r=l,2 ^ j=-l 
' r=l,2 ^ j=-l 



(A.7) 



where 

Ij = & dzexp(-(3z)Gjj(z), (A. 8) 

and for any couple of integers or zeroes (J, k), 

Gj : k(z) =< n + j\{z - /jO^K + k> (A.9) 

is the Green function connecting the eigenstates |n +j > 
and \no + k > of h^. The index r appearing in Eq. i|A.7|) 
is the sub-lattice index. The operator (z — hi)~ l appear- 
ing in the right-hand side of Eq. ljA.91) is the so-called 
"resolvent operator" . To evaluate the Green functions 



Gj t k(z), we have to know how the resolvent operator acts 
on the ket \n + k >. According to Eq. (|A.3|I . the action 
of the Hamiltonian operator hi is nothing but the action 
of the operator h^ plus the action of the operator hi. 
Their action can be determined explicitly as follows. If 
A and B are two operators, the following identity holds: 

i - i = = . (a.io) 

ABA B B A 

Then, with the identifications A = z — hi and B = z — h^, 
one has 



z — hi 



- hi 



-hj- 



-hi 



-h L 



-hj- 



- hL z- hi 



(All) 
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The right-hand side of Eq. (jA.llfl will be useful in the 
evaluation of the propagators Gj k{z) in the basis formed 

by the eigenstates of Hl- In principle, for each values of 
j , one needs to construct a p x p system of equations in 
the variables Gj i / C (z), where the order p of each system 
is equal to the cardinality of the chosen basis. On the 
other hand, only those functions Gj^(z) connecting ba- 
sis vectors will give non-vanishing contributions. Hence, 



in our case we have to deal with four systems, each of 
these made up of four equations. We will write down the 
explicit form of such a system for the variables Gjj(z) 
that, as we can see from the Eq. (|A.7|) . are the needed 
quantities to determine the free energy per site. In each 
sub-lattice labeled by r, and omitting the index for the 
function Gjj(z) , we have 



6j,o = (z + (S-K X r))G ,o(z) + 



k VK + iKW+2) 



G2,2( Z ) VV 



Sj,i = (z-(8-K X r))G 1A (z) + 



G ,o( z ) 



J F > A K^/noino + l) ri 
H 2 G 2,2(z) <j> 2 G-i-i(z) ip r ; 

(z-Uo + 3(5-K X r))G-i,-i(z) 



7T G ,o{z) (p G 1A {z) 1p r 



5j,2 - {z-U -3(5-Kxr))G 2 ,2{z) + 
Ky/(no + l)(no + 2) 

7. Go,o(Z) Ipr 



Gi,i(z)<f> 



(A.12) 



r 



Therefore, when j = 0, j = 1, j — —1, and j = 2 
solving system Eq. 1|A.12() provides, respectively, Go,o{z), 
Gx t i(z), G_i._i(z), and G 2 ^(z) in each sub-lattice. Each 
of these solutions may be written as 



Gj, 3 (z) 



Nj(z) 
D(z) 



(A.13) 



Since we are operating in the strong-coupling regime, 
we retain only the contributions proportional to non- 
vanishing powers of Uq in Nj(z) and D(z). We can now 
obtain the explicit expression for the free energy by direct 
evaluation of the second line of Eq. (|A.7(1 . First, the in- 
tegrals Ij appearing in Eq. (|A.8(1 are solved by the usual 
integration techniques in the complex plane, and we de- 
termine the poles of the functions Gjj(z) by solving the 
equation 



The roots of Eq. I|A.14|) provide the eigenvalues of the 
Hamiltonian hi at the needed order of approximation 
and allows to calculate explicitly the right-hand side of 
Eq. (|A.8p . Finally, inserting the expressions for the quan- 
tities Ij in the second line of Eq. (|A.7|I . yields the desired 
expression for the free energy per site Eq. (|3U|I . 

The method of the resolvent allows as well to obtain 
the explicit expression for the mean field order parame- 
ter x m eacn sub-lattice. Following the same procedure 
adopted to evaluate the free energy per site, the "mean 
number" x m a given sub-lattice can be determined by 
the the formula 



D(z) = 



(A.14) 



X 2 



En < n\hexp(-f3hj)\n > _ + 1 ' 
J2„ < n\exp(-0hi)\n > ° 2 - 
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no/o + (no + l)/i + (no-l)/-i + (n + 2)7 2 1 
I +I 1+ I_ 1+ I 2 ( "° + 2 } 



(A.15) 



r 



where the integer index n runs over the finite number of the generic sub-lattice r finally reads 
local eigenstates being considered. The magnetization of 



J 



Xr ~ 2 



If 5-KXr , r n(x « (J F < A <P) 2 (n + l) 

2\~ r tanh ^ K + TT ^-- 



-(J F - A <j>) 2 ) tanh[/3(A r + ^)] 



where 



2U Xf 



2K (n + 1)(5 - K Xr ) 4? $ + K 2 ijj 2 (n + 1) 



(A.16) 



a r 



2(J^ A (/ ) yKMno + ir + ((J 



F,A±\2 t^2„/,2 



r 



K z ri(n + l))(5-K Xr ) 



(A.17) 



In the ferromagnetic-like case, the single atom hopping 
amplitude is J F , so that \\ = \2 = X, and the two 
subdattices are characterized by the same magnetiza- 
tion, that is, by the same constant filling factor n . In 
the antiferromagnetic-like case, the single atom hopping 



amplitude is J A , so that xi = ~ Xi = Xi an d the two 
sub-lattices have opposite magnetizations, that is, two 
different constant filling factors no and n + 1, respec- 
tively. 
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